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ABSTRACT 

We implement the hydro-PM (HPM) technique (Gnedin & Hui 1998) in the 
hydrodynamical simulation code GADGET- 2 and quantify the differences be- 
tween this approximate method and full hydrodynamical simulations of the 
Lyman-Q forest in a concordance ACDM model. At redshifts z — 3 and z = 4, 
the differences between the gas and dark matter (DM) distributions, as mea- 
sured by the one-point distribution of density fluctuations, the density power 
spectrum and the flux power spectrum, systematically decrease with increas- 
ing resolution of the HPM simulation. However, reducing these differences to 
less than a few percent requires a significantly larger number of grid-cells than 
particles, with a correspondingly larger demand for memory. Significant differ- 
ences in the flux decrement distribution remain even for very high resolution 
hydro-PM simulations, particularly at low redshift. At z = 2, the differences 
between the flux power spectra obtained from HPM simulations and full hy- 
drodynamical simulations are generally large and of the order of 20-30%, and 
do not decrease with increasing resolution of the HPM simulation. This is due 
to the presence of large amounts of shock-heated gas, a situation which is not 
adequately modelled by the HPM approximation. We confirm the results of 
Gnedin & Hui (1998) that the statistical properties of the flux distribution 
are discrepant by <; 5 — 20% when compared to full hydrodynamical simula- 
tions. The discrepancies in the flux power spectrum are strongly scale- and 
redshift-dependent and extend to large scales. Considerable caution is needed 
in attempts to use calibrated HPM simulations for quantitative predictions of 
the flux power spectrum and other statistical properties of the Lyman-a forest. 

Key words: Cosmology: intergalactic medium - large-scale structure of uni- 
verse - quasars: absorption lines - hydrodynamics - methods: numerical 



1 INTRODUCTION 

The prominent absorption features blue-ward of the 
Lyman- a emission in the spectra of high-redshift 
quasars (QSOs) are believed to arise from smooth den- 
sity fluctuations of a photoionised warm intergalactic 
medium which trace the dark matter distribution in a 
relatively simple manner (see Rauch 1998 and Wein- 
berg et al. 1999 for reviews). As a result, the flux power 
spectrum of this 'Lyman-a forest' has become a pow- 
erful quantitative probe of the matter power spectrum 
on scales of 1 h' 1 to 40 h" 1 Mpc at redshifts z = 2 — 4. 
At these scales and redshifts, the matter distribution is 
linear or mildly non-linear, a regime that can be accu- 
rately modelled with numerical simulations. Such simu- 



lations have been used to obtain quantitative estimates 
of the clustering amplitude and constraints on cosmo- 
logical parameters from the Lyman-Q forest (Croft et 
al. 1998; Croft et al. 1999; McDonald et al. 2000; Hui 
et al. 2001; Croft et al. 2002; McDonald 2003; Viel et 
al. 2003; Viel, Haehnelt & Springel 2004; Viel, Weller 
& Haehnelt 2004; McDonald et al. 2004a; Desjacques & 
Nusser 2004) or on astro-physical parameters (Theuns 
et al. 1998, Meiksin et al. 2001, McDonald et al. 2004b, 
Bolton et al. 2005). 

Unfortunately, the flux power spectrum does not 
only depend on the dark matter (DM) distribution but 
also on the thermal state of the intergalactic medium 
(IGM), and possibly on feedback effects due to star for- 
mation and active galactic nuclei (AGN). Ideally, one 
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would like to use simulations which not only take into 
account the non-linear gravitational clustering of the 
matter distribution but also all the relevant hydrody- 
namics of the gas, including effects of galaxy forma- 
tion physics, such as radiative cooling and heating, star 
formation and winds driven by stellar associations or 
AGN. However, full hydrodynamical simulations of the 
Lyman-a forest are computationally very demanding. 
This makes their use for extensive parameter studies 
difficult. In addition, some physical processes, such as 
the feedback mechanisms, are still poorly understood. 
Thus, the use of approximate numerical calculations of 
the flux distribution of the Lyman-a forest very attrac- 
tive, an approach that has been widely applied in previ- 
ous work (e.g. McGill 1990; Hui et al. 1997; Meiksin & 
White 2001; Viel et al. 2002b; Zhan et al. 2005). Note 
that such approximate calculations of the Lyman-a flux 
distribution have been crucial in establishing the mod- 
ern paradigm for the origin of the Lyman-a forest in 
the first place (Bi 1993, Bi & Davidsen 1997, Viel et al. 
2002a). 

In 1998, Gnedin & Hui (GH) have proposed the 
'Hydrodynamic Particle-Mesh method' (HPM) as an ef- 
ficient numerical method to approximate the formation 
and evolution of the Lyman-a forest. This technique 
is based on a particle-mesh (PM) approach for follow- 
ing the evolution of dark matter. The gravitational po- 
tential of the PM solver is then modified with an ef- 
fective potential which mimics the effect of gas pres- 
sure. GH found that global statistical properties of the 
flux distribution in HPM simulations are accurate to 
~ 5 — 20% when compared to full hydrodynamical sim- 
ulations. This prompted e.g. McDonald et al. (2004a) to 
use HPM simulations that were calibrated with a small 
number of hydrodynamical simulations to obtain pre- 
dictions of the flux power spectrum for a wide range 
of cosmological and physical parameters describing the 
thermal state of the gas. 

The statistical errors of the flux power spectrum 
obtained from high-resolution Echelle spectra are ~ 4 % 
and can in principle become as small as a few percent 
for large samples of low-resolution spectra (e.g. Kim et 
al. 2004, McDonald et al. 2005). This has opened up the 
exciting prospect to use the Lyman-a forest to constrain 
inflationary parameters and the nature of dark matter, 
based on high accuracy measurements of the DM power 
spectrum inferred from the Lyman-a forest (Viel, Weller 
& Hahenelt 2004; Seljak et al. 2004; Viel et al. 2005). 
However, a prerequisite is the availability of accurate 
predictions of the flux power spectrum for a wide range 
of parameters. 

The hydrodynamical code GADGET- 2 (Springel, 
Yoshida & White 2001; Springel 2005), which we have 
used extensively in earlier work for full hydrodynami- 
cal simulations of the Lyman-a forest (Viel, Haehnelt 
& Springel 2004; Bolton et al. 2005), is a TreeSPH code 
which also offers a PM algorithm which can optionally 
be used to calculate long-range gravitational forces. In 
this code the HPM method of GH can therefore be easily 
implemented. This makes GADGET-2 well suited for a 
detailed analysis of the accuracy and systematic uncer- 
tainties of the HPM method by comparing simulations 



run with it to full hydrodynamical TreeSPH-PM simu- 
lations. 

In this paper, we perform such an analysis and in- 
vestigate the dependence of the discrepancies between 
HPM and full hydrodynamical simulations on a range 
of numerical parameters for the relevant redshift range 
2 = 2 — 4. Note that we here do not intend to optimize 
the HPM method. 

The outline of the paper is as follows. In Sec- 
tion 2, we briefly describe the hydrodynamical code 
GADGET-2 and we review the basic equations of the 
HPM formalism. We also show that the HPM imple- 
mentation of GADGET-2 and the HPM code by Gnedin 
& Hui give similar results for a suitable choice of numer- 
ical parameters. In Section 3, we discuss the differences 
between our HPM implementation and full hydrody- 
namical simulation by analysing the statistical proper- 
ties of the flux distribution. We further analyse the ef- 
fect of shock heating, the influence of various numerical 
parameters on the results, and the CPU time and mem- 
ory requirements. Finally, Section 4 contains a summary 
and our conclusions. 



2 SIMULATION METHODS OF THE 
Lyman-a FOREST 

2.1 Full hydrodynamical simulations 

The hydrodynamical simulation code GADGET-2 
(Springel, Yoshida & White 2001; Springel 2005) can op- 
tionally employ a PM technique to calculate long-range 
gravitational forces, resulting in a 'TreePM' scheme for 
gravitational forces. We will use hydrodynamical sim- 
ulations run with this SPH/TreePM implementation of 
GADGET-2 as "reference" simulations to assess in detail 
the accuracy and systematic uncertainties of the approx- 
imate HPM method. The TreePM approach speeds up 
the calculation of long-range gravitational forces consid- 
erably compared to a tree-only implementation. 

All our simulations were performed with periodic 
boundary conditions and an equal number of dark 
matter and gas particles. We employ the 'entropy- 
formulation' of SPH proposed by Springel & Hernquist 
(2002). Radiative cooling and heating processes are fol- 
lowed using an implementation similar to that of Katz 
et al. (1996) for a primordial mix of hydrogen and he- 
lium. We have assumed a mean UV background pro- 
duced by quasars as given by Haardt & Madau (1996), 
which leads to reionisation of the Universe at z ~ 6. 
The simulations are run with heating rates increased by 
a factor of 3.3 in order to achieve temperatures which 
are close to observed temperatures (Abel & Haehnelt 
1999, Schaye et al. 2000, Ricotti et al. 2000). 

In order to maximise the speed of the dissipative 
hydrodynamical simulations we have employed a sim- 
plified star-formation criterion in the majority of our 
runs. All gas at densities larger than 1000 times the 
mean density was turned into collisionless stars. The 
absorption systems producing the Lyman-a forest have 
small overdensity so this criterion has little effect on flux 
statistics, while speeding up the calculation by a factor 
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of ~ 6, because the small dynamical times that would 
otherwise arise in the highly overdense gas need not to 
be followed. In a pixel-to-pixel comparison with a sim- 
ulation which adopted the full multi-phase star forma- 
tion model of Springel & Hernquist (2003) we explicitly 
checked for any differences introduced by this approxi- 
mation. We found that the differences in the flux proba- 
bility distribution function were smaller than 2%, while 
the differences in the flux-power spectrum were smaller 
than 0.2 %. We have also turned off all feedback options 
of GADGET-2 in our simulations. An extensive resolu- 
tion and box size study has been performed in Viel, 
Haehnelt & Springel (2004) and in Bolton et al. (2005). 

For all simulations presented here we have adopted 
a box size of 30 comoving h~ Mpc and the cosmolog- 
ical parameters f2om = 0.26, Hoa = 0.74, Q oh = 0.0463 
and H = 72 km s _1 Mpc _1 , cr 8 = 0.85 and n = 0.95 
(the parameters of the B2 simulation in Viel, Haehnelt 
& Springel 2004). The CDM transfer functions of all 
models have been taken from Eisenstein & Hu (1999). 

2.2 HPM implementation of GADGET-2 

GH proposed to introduce an effective potential that 
mimics gas pressure into an otherwise collisionless 
dark matter simulation, carried out with a particle 
mesh code. This method has become known as Hydro- 
Particle-Mesh (HPM) approximation. The idea of the 
HPM approximation is to take advantage of the fact 
that the low density IGM responsible for most of the 
Lyman-a forest absorption obeys a simple relation be- 
tween gas density and gas temperature, which is well 
described by a power-law 'equation of state': 

T = T (z) (1 + J)^" 1 . (1) 

The evolution of To and 7 with redshift depends on the 
reionisation history (Hui & Gnedin 1997). The 'equation 
of state' predicts the temperature of gas of given density 
to better than 10% for the low density IGM where shock 
heating is not important. Instead, the temperature is 
set by a balance between photoionisation heating and 
adiabatic cooling due to the expansion of the universe. 

Based on the density alone, equation (1) also allows 
an estimate of the thermal pressure which enters the 
equation of motion for a cosmic gas element. We know 
from full hydrodynamical simulations that the baryons 
follow the dark matter generally well apart from high 
density regions where pressure effects on small scales 
become important. GH suggested therefore to use the 
density of the dark matter in a PM simulation together 
with Eqn. (1) to estimate the temperature and pressure 
of the gas. One can then obtain the acceleration on a 
cosmic gas element due to the gradient of the pressure 
as 

*L + Hv = -V4>-±VP, (2) 

where v is the gas peculiar velocity, 4> is the gravitational 
potential, and P is the thermal pressure. If the gas is 
highly ionised (so that the mean molecular weight is 
roughly constant, which is true for the Lyman-alpha 
forest), and the temperature is a function of density 



3 

only, so that P — P(p), equation (2) can be reduced to 
the expression 

!^ + Hv = -V^, (3) 
where 

i> = <t> + n, (4) 

and H, the specific enthalpy, is 

Equation (3) is identical to the equation of motion 
for the collisionless dark matter except that the usual 
gravitational potential <f> is replaced by an effective po- 
tential if), which takes into account both gravity and 
thermal pressure. Since the gravitational potential (f> has 
to be computed from the density field in a regular PM 
simulation anyway, computing the enthalpy adds only a 
modest computational overhead. 

We have implemented this HPM method in the sim- 
ulation code GADGET-2. We closely follow the approach 
of GH with only a few minor differences. In the HPM 
code of GH, only one set of particles was used, i.e. the 
fact that the dark matter does not feel the pressure on 
small scales was neglected. As GADGET-2 is a SPH code 
which treats DM and baryons separately, we kept this 
distinction in our HPM implementation. This may re- 
sult in some small differences on small scales. In Section 
2.3, we will compare simulations with the HPM imple- 
mentation of GADGET-2 to runs carried out with the 
HPM code of GH (kindly provided by Nick Gnedin). 

There are three numerical parameters defining 
the technical details of our HPM implementation in 
GADGET-2 . The first parameter is the number of cells 
of the PM grid. We describe this by iVgrid, the num- 
ber of cells per dimension. The second parameter, H s , 
describes the scale of the smoothing applied to the en- 
thalpy field before taking its spatial derivative. The den- 
sity and enthalpy fields are more sensitive to shot noise 
than the gravitational potential, because for the latter, 
high frequency noise is suppressed as <f>(k) oc 5 k k~ 2 . We 
have thus followed GH and apply a Gaussian smooth- 
ing to the density field before computing the enthalpy 
and its spatial derivative. We apply a smoothing fac- 
tor exp(— k 2 h 2 ) to the density field in Fourier space, 
where h s = H s L/N gI id- The third numerical parameter, 
r s = A s L/N gTi d, is the scale of the smoothing of the PM 
force, which we usually express in terms of A B , i.e. in 
units of the mesh cell size. The parameter ^4 S hence 
controls the level of residual force anisotropies in the 
PM force. In the TreePM code, r B also gives the scale of 
the short-range/long-range force split. We will discuss 
the choice of numerical values for these parameters in 
Section 3.3. Note that the HPM code of GH has only 
two parameters Af gr id and H B , i.e. no attempt is made 
to make the PM force more isotropic on the scale of the 
mesh. GH have adopted the choice H B — 3. 

To fix the slope and normalisation of the power- 
law temperature-density relation of the IGM, our code 
follows the thermal history of two fiducial gas elements 
at density values equal to the mean cosmic density and 
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Figure 1. Power spectrum of the gas density field of the HPM simulations run with GADGET-2 at z = 3, at two different 
resolutions and for several different values of the parameter iVg r jd- The power spectrum of the full hydrodynamical simulation 
is represented by the filled triangles. 
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Figure 2. Left: differences in the probability distribution functions of the dark matter density field between GADGET-2 (G2) 
and the Gnedin & Hui (GH) code. Both of them have been run in the PM mode with a grid of 200 3 for GH and 400 3 for G2. 
Right: Fractional differences in the 3D matter power spectrum. The results are shown at three different redshifts z = 2, 3, 4 as 
dashed, continuous and dotted curves, respectively. 



at 1.1 times the mean cosmic density. For a specified 
evolution of the ionising UV background, we can then 
compute the values of To and 7 from the temperatures 
attained by these two fiducial gas elements. 

In Figure 1, we compare the 3D gas power spectrum 
for a range of HPM simulations with different parti- 
cle numbers and mesh sizes with a full hydrodynamical 
simulation with 200 3 dark matter and 200 3 gas parti- 



cles (shown as triangles). All simulations were run with 
GADGET-2 . We only show results at z = 3, but note 
that the results at z = 2 and z — 4 are very similar. On 
large scales (k < 6/1/Mpc), the power spectrum of the 
gas distribution of HPM simulations converges nicely 
to that of the full hydrodynamical simulation when the 
resolution of the mesh used for calculating the gravi- 
tational forces is increased. Note, however, that even 
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Figure 3. Left: differences in the probability distribution functions of the gas density field between simulations run with 
GADGET-2 (G2) and the Gncdin & Hui (GH) code. Both of them have been run in the HPM mode with a grid of 200 3 for 
GH and 400 3 for G2. Right: Fractional differences in the 3D matter power spectrum. The results are shown at three different 
rcdshifts z = 2, 3, 4 as dashed, continuous and dotted curves, respectively. 
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Figure 4. Left: differences in the probability distribution functions of the dark matter density field between simulations run 
with GADGET-2 in its HPM and in its TreePM mode (the PM grid for the TrcePM run is fixed to the value iV grid = 200). 
Right: Fractional differences in the 3D matter power spectrum. The results are shown at z = 3 and for three different values of 
iVgrid (400,600,1200) as continuous, dashed and dot-dashed curves, respectively. 



for very high resolution (six times more mesh cells in 
the HPM simulation than particles in the full hydro- 
dynamical simulation) the power on small scales in the 
HPM simulations is significantly smaller than that in 
the full hydrodynamical simulations. Note further that 
changing the mesh resolution is more important than 
changing the particle number in the HPM simulations. 
The thin and thick solid curves are for HPM simula- 
tions with the same grid resolution but a factor eight 



different particle number. They are virtually identical. 
We also note that the results and trends for the dark 
matter power spectrum are qualitatively similar. In the 
runs discussed in the following, we will use the HPM 
implementation of GADGET-2 with 2 x 200 3 particles 
and with Af gri d > 200. 
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Figure 5. Left: differences in the probability distribution functions of the gas density field between simulations run with 
GADGET-2 in its HPM and in its TrccPM mode (the PM grid for the TreePM run is fixed to the value 200). Right: Fractional 
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as continuous, dashed and dot-dashed curves, respectively. 



2.3 Comparison between the HPM 

implementation of GADGET-2 and the 
HPM code of Gnedin & Hui 

In this section we compare the gas and dark matter dis- 
tribution of simulations run with the HPM implementa- 
tion of the GADGET-2 code and the HPM code of GH. 
We use the same initial conditions and temperature- 
density relation. At z = 2, 3, and 4, To and 7 (To, 7) have 
the following values: (21500 K,1.505), (21500 K,1.524) 
and (19200 K,1.536). 

In Figures 2 and 3, we show the relative differ- 
ences of the probability distribution and the power spec- 
trum of the dark matter and gas density at redshifts 
2 = 2, z — 3, and 2 = 4. We have varied the reso- 
lution of the mesh to calculate the gravitational force 
in the HPM implementation of GADGET-2 in steps of 
factors of two. The other two relevant parameters in 
the GADGET-2 runs have been set to H s = 3 and 
A s — 1.25. For the case shown in Figures 2 and 3, 
the grid resolution for the HPM implementation of 
GADGET-2 was a factor two higher than that used for 
the HPM code of GH. In this case the agreement was 
best, better than 5% (dark matter) and 8% (gas) for 
the probability distribution function* (pdf) and better 
than 2% for power spectra at wavenumbers relevant for 
constraining cosmological parameters with the Lyman- 
a forest, 0.3 < k (ft/Mpc)< 3 (Viel, Haehnelt & Springel 
2004). Because of the smoothing applied to the PM 
force in GADGET-2 , a somewhat finer mesh is needed 
to match the results of the HPM code by GH, where 



The pdf is defined as the number of points or pixels in a 
given x-axis bin with the property that its integral along the 
x-coordinate is one. 



(CODE, # part.) 


Afgrid 


CPU-time (ks) 


Mem. (Gb) 


HPM-200 3 


100 


1.5 


3 


HPM-200 3 


200 


3.1 


3.5 


HPM-200 3 


400 


4.7 


4.5 


HPM-200 3 


600 


11.2 


12 


HPM-200 3 


1200 


15 


76 


HPM-400 3 


100 


33 


26 


HPM-400 3 


200 


35 


28 


HPM-400 3 


400 


40 


30 


HPM-400 3 


600 


44 


36 


Hydro-200 3 


200 


183 


3.2 


Hydro-400 3 


400 


11700 


26 


GH HPM-200 3 


200 


5.4 


3.5 



Table 1. CPU-time required to reach 2 = 2 for simulations 
of a 30 Mpc/h box ACDM model and for several different 
resolutions and values of the parameter N SI ^. The memory 
required is shown in the last column. All the values are wall- 
clock times for 32 CPUs (1.3 GHz Itanium 2) of the SGI 
Altix 3700 (COSMOS) at DAMTP (Cambridge). 



such a smoothing is not carried out and larger force 
anisotropies on the mesh scale are accepted. By reduc- 
ing A B , the agreement of the two codes could be im- 
proved further. The two HPM codes agree very well. In 
the following we will only use the HPM implementation 
of GADGET-2 but our results should apply similarly to 
the GH code. 
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3 COMPARISON BETWEEN FULL 
HYDRO DYNAMICAL AND HPM 
SIMULATIONS 

3.1 The dark matter and gas density fields 

We first want to check the agreement of the dark mat- 
ter and gas distributions between simulations run with 
the TreePM and HPM implementations of GADGET-2 . 
In Figure 4, we show the differences in the density pdf 
and the power spectrum for the dark matter distribu- 
tion at z — 3, for three different values of Afgrid used in 
the HPM simulation. The results at z — 2 and z — 4 
are similar. The simulations were run with 200 3 and 
2 x 200 3 particles, respectively. In the simulation with 
the Tree-PM implementation, the number of mesh cells 
of the PM grid was set equal to the number of particles. 
As also expected from the results shown in Fig. 1, the 
differences become smaller with increasing resolution of 
the PM grid used for the HPM implementation. The dif- 
ferences in the pdf of the DM density are smaller than 
10% (20%) for AT grid =600 (400). If a very fine mesh of 
dimension A r gr id = 1200 is used, the pdf of the HPM sim- 
ulation is indistinguishable from that of the full hydro- 
dynamical TreePM simulation. For AT grid =600 (400) the 
discrepancy in the dark matter power spectrum (right 
panel) is less than 2% (4%) for 0.2 < fc(/i/Mpc) < 2. For 
A?g r id=1200, the difference is less than 0.5% in the same 
range of wavenumbers. At larger wavenumber the differ- 
ences in the power spectra become much larger due to 
the much higher resolution achieved with the TreePM 
code. Note, however, that these small scales are not used 
for the recovery of the dark matter power spectrum from 
the Lyman-a forest because of the uncertainties in the 
flux power spectrum due to the thermal history and the 
metal contamination of the IGM (Kim et al. 2004) . 

Figure 5 shows the difference in the gas distribu- 
tions between simulations with the HPM and TreePM 
implementations. The differences are similar to those 
found in the dark matter distribution. 



3.2 Flux statistics 

3.2.1 The flux probability distribution function 

The flux distribution in the Lyman-a forest depends on 
the spatial distribution, the peculiar velocity field and 
the thermal properties of the gas. In the last section, we 
have shown that the gas distribution of the HPM sim- 
ulations converges rather well to that of the full hydro- 
dynamical simulations when the resolution of the PM 
mesh is improved. For the flux distribution the situa- 
tion is more complicated, however. In Figure 6, we plot 
the differences in the pdf of the flux for HPM simula- 
tions with a range of A gr id values compared with the 
full hydrodynamical simulations at z = 2, 3, 4. The sim- 
ulations are the same as those discussed in section 3 and 
shown in figures 4 and 5 (these figures show results only 
at z = 3). The curves without symbols show the results 
for the same amplitude of the ionising UV background 
as in the full hydrodynamical simulations. Note that 
this means that the the flux distribution has not been 
rescaled to a fixed mean flux, as it is often done. Such 



a rescaling would mask the numerical effects we seek 
to identify here. However, to facilitate comparison with 
other work (e.g. McDonald et al. 2004a), the curves 
with triangles show the pdf of the flux after re-scaling 
the flux distribution of the N slid = 1200 HPM simula- 
tion such that the mean transmitted flux is the same as 
in the full hydrodynamical simulations. 

At z = 3, the flux distribution of the HPM simu- 
lations converges reasonably well to that of the full hy- 
drodynamical simulations. With the exception of flux 
levels F > 0.8, the differences are smaller than 5% for 
A^rid = 600 and even smaller for higher resolutions of 
the PM mesh. In regions of low absorption (F > 0.8) 
the differences are, however, large (10-20%), change sign 
with increasing resolution, and do not converge. We 
have inspected a few spectra individually and found that 
the discrepancy is due to differences in both density and 
temperature in the lowest density regions. At z — 4 
these differences in regions of low absorption are sub- 
stantially larger. Because of the strong decrease of the 
mean flux with increasing redshift, these regions corre- 
spond to significantly more underdense regions than at 
z = 3. At z — 2 additional large differences up to 50% 
arise in regions of strong absorption, which also do not 
vanish with increasing resolution. For the iV gr id = 1200 
HPM simulation, the overall agreement with the full 
hydrodynamical simulation is of the order of 2% for 
F < 0.85 at z = 3, 4, while at z = 2, discrepancies 
of the order of > 10% remain both in underdense and 
very dense regions. The differences at z — 2 and for 
F < 0.15 are due to the gas in dense regions being sub- 
stantially colder in the HPM simulations than in the full 
hydrodynamical simulation where a significant portion 
of the dense gas is shock heated. In Figure 6 we over- 
plot the results from a higher resolution HPM run with 
2 x 400 3 particles and A" grid = 600 as a long dashed line 
with filled diamonds. The results are very similar to the 
HPM simulation with 2x200 3 particles and AT grid = 200. 
We hence confirm the findings of GH that the differences 
in the flux pdf between HPM and full hydrodynamical 
simulations are of the order of 10-15%. 

3. 2. 2 The flux power spectrum 

The main motivation of the use of HPM simulations 
comes presently from the need for accurate predictions 
of the flux power spectrum for a wide range of astrophys- 
ical and cosmological parameters. Such a grid of predic- 
tions allows a detailed comparison with observational 
data and a determination of best-fit values and confi- 
dence intervals of cosmological parameters (McDonald 
et al. 2004a). 

In Figure 7, we plot the differences of the flux power 
spectrum of HPM simulations with a range of mesh 
sizes compared with full hydrodynamical simulations 
at 2 = 2, 3, 4. The simulations are the same as those 
discussed in the previous sections. As in figure 6, the 
curves without symbols show results for the same am- 
plitude of the ionising UV background while the curves 
with empty triangles show the flux power spectrum af- 
ter rescaling the flux distribution of the Ar grid = 1200 
HPM simulation such that the mean flux is the same 
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Figure 6. Effect of the parameter N gI [^ for simulations of a 30 Mpc/h box with 2 X 200 (gas and dm) particles. Left 
panel: Fractional differences between the probability distribution functions of simulations with A r gr ; ( j=200 (dotted), Wg r i<j=400 
(continuous), JV gr ;d=600 (dashed), and -/V gr ;d=1200 (dot-dashed) at z = 2. Also shown is the full hydrodynamical TreePM 
simulation with 7V gr id=200 and with the same initial conditions. The long-dashed line with filled diamonds represents results 
for the higher resoultion run with 2 X 400 3 particles and with 7V gr ; t j=600. The other HPM parameters have been fixed to the 
fiducial values described in the text. The dot-dashed curve with ovcrplottcd empty triangles is for a simulation with A r gr i ( j=1200 
for which the simulated flux has been scaled to match the value of the full hydro simulations, in the other cases the spectra 
have not been scaled (see text for the details). Middle Panel: Results at z = 3. Right Panel: Results at z = 4. 



as in the full hydrodynamical simulations. In Figure 7 
we show the results from a higher resolution HPM run 
with 2 x 400 3 particles and N gI id = 600, as the long- 
dashed line with overplotted filled diamonds. At red- 
shift z = 4 and 2 = 3 there is perfect agreement with 
the Af gr id = 1200 HPM simulation, in the wave number 
range of interest here. At z — 2 there are small dif- 
ferences of the order of < 5%. Thereby, increasing the 
number of particles does not improve the agreement sig- 
nificantly. 

At redshifts z = 3 and z = 4, the flux power spectra 
of the HPM simulations converge well to those of the full 
hydrodynamical simulations, but only for resolutions of 
the PM mesh where the number of the mesh cells is sub- 
stantially larger than that of the number of particles in 
the full hydrodynamical simulations. At z — 3, the HPM 
simulations with AT gr id=400 (600) have scale- dependent 
differences of about 10% (7%) in the wavenumber range 
relevant for inferring the matter power spectrum. For 
A r gr id=1200, there is a scale-independent offset of about 
5% (3% when rescaled to the same mean flux). At red- 
shift z — 4 the situation is very similar. However, at 
redshift z = 2, the flux power spectrum of the HPM 
simulations does not converge to that of the full hydro- 
dynamical simulation. The differences are here actually 
smallest for the HPM simulation with lowest resolution 
(A r gr id=200). However, even in this case the discrepan- 
cies are large and strongly scale dependent, of the order 
of 25-30% at the largest scales. At small scales k > 0.02 
s/km, the size of the disagreement and its scale de- 
pendence is similar to that found by McDonald et al. 
(2004a, their figure 5). Note that because of the smaller 



box size of their hydro simulations, McDonald et al. were 
not able to probe scales k < 0.007 s/km (at z — 2), 
where the differences increase dramatically. Note that 
the amount of shock-heated gas is significantly larger in 
simulations with larger box-size. To test further to what 
extent these discrepancies at large scales depend on the 
resolution of the hydro-simulation, we have run an addi- 
tional hydro-simulation with 64 times higher mass reso- 
lution (2 x 200 3 particles in a 7.5 h/ Mpc box). There is 
good agreement with the results shown in Figure 7. We 
stress here that our goal is to get a good convergence of 
the flux power in the range 0.003 < k (s/km) < 0.03, 
which is the range which is used for the matter power 
spectrum reconstruction as in Viel et al. (2004). 

These large differences and the lack of convergence 
appear perhaps counterintuitive considering the rather 
good convergence of the gas and dark matter distribu- 
tion. However, they simply originate in the large differ- 
ences in the pdf of the flux distribution, which in turn 
are due to the different thermal state of the gas in high 
density regions in the HPM and full hydrodynamical 
simulations. 

At redshift z = 2, a larger proportion of the ab- 
sorption is from gas in high density regions, which is 
shock-heated in the full hydrodynamical simulations 
and therefore on average hotter than in the HPM sim- 
ulations. This tends to mainly affect the strong absorp- 
tion systems which contribute significantly to the flux 
power spectrum at large scales (Viel et al. 2004a; Mc- 
Donald et al. 2004b). We fill discuss this further in the 
next section. 
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Figure 7. Effect of the parameter N SI ^ for simulations of a 30 Mpc//i box with 2 X 200 3 (gas and dm) particles. Left panel: 
Fractional differences between the ID flux power spectra of simulations with A r gr ; ( j=200 (dotted), A r gr i < j=400 (continuous), 
A r gr ; ( j=600 (dashed), and A r gr ; t i=1200 (dot-dashed) at z = 2. Also shown is the full hydrodynamical TreePM simulation with 
Af gr i(j=200 and the same initial conditions. The other HPM parameters have been fixed to the fiducial values described in the 
text. The dot-dashed curve with overplotted empty triangles is for a simulation with Af gr ; t i=1200 for which the simulated flux 
has been scaled to match the value of the full hydro simulations, in the other cases the spectra have not been scaled (see text for 
the details). The long-dashed line with filled diamonds represents results for the higher resoultion run with 2 X 400 3 particles 
and with Af gr ; t j=600 (results scaled to reproduce the same T c ff). Middle Panel: Results at z = 3. Right Panel: Results at z = 4. 
In all the panels the dashed area represents the range of wavenumbers used by Viel, Haehnelt &; Springel (2004) to recover 
cosmological parameters. 



3.2.3 Temperature effects on the flux pdf and the flux 
power spectrum 

We have argued that the approximation of the relation 
between gas density and gas temperature as a power- 
law breaks down at low redshift. This approximation 
inevitably does not take into account the amount of 
moderately shock-heated gas that is falling into the po- 
tential wells of the dark matter haloes. In this section, 
we want to check this explicitly. For this purpose we use 
the hydrodynamical simulation and the HPM run with 
jV grid =600. 

As a first step we perform the following simple test. 
We superimpose onto the full hydro-dynamical SPH 
simulation the temperature-density relation of the HPM 
runs, and then recompute the QSO spectra. We find 
that this results in differences much smaller than those 
in Figure 7, of order 8%, 5% and 4%, at z = 2, 3 and 
4, respectively, at the largest scales. Most of the dis- 
crepancy is thus indeed due to the differences in the 
thermal state especially at low redshift. Differences in 
the thermal state will lead, however, also to pressure 
differences during the dynamical evolution, which will 
modify the mass distribution and the peculiar velocity 
field. In shock fronts, the change in particle trajecto- 
ries can be substantial. Since the HPM implementation 
does not capture shocks, it would not treat the dynam- 
ics correctly even if the temperatures would be accurate 
at all times. To investigate this further we have run an 
SPH simulation with artificial viscosity set to zero and 



the temperature-density relation of the HPM simula- 
tion. This should mimick an 'ideal' HPM simulation: 
the gravitational force is resolved with high accuracy 
and in an isotropic way, while the pressure gradients 
are smooth and resolved everywhere with the maximum 
resolution allowed by the local particle sampling. The 
standard HPM method has a less well resolved gravita- 
tional force and should be sensitive to over- or under- 
smoothing of the pressure field in regions of high or 
low particle density, respectively. The results are shown 
in Figure 8. The SPH simulation is represented by the 
dashed line while the dotted line is for the HPM simu- 
lation with iVg r id — 1200 (both the runs have the same 
number of particles equal to 2 x 200 3 ). There is good 
agreement with the HPM simulations, suggesting that 
the discrepancy in the flux power is primarily due to 
the different thermal state of the gas due to shocks and 
not to any artefacts of our particular HPM implementa- 
tion. The total effect on the flux power spectrum should 
thereby be a combination of an increase of the over- 
all amount of shock-heated gas with decreasing redshift 
and the change of the mean effective optical depth. The 
flux power spectrum becomes increasingly sensitive to 
higher-density gas with decreasing redshift due to the 
decreasing effective optical depth. 

We will now investigate the relation between the 
differences between HPM and SPH and gradients in the 
velocity field of the gas. Negative gradients in the pecu- 
liar velocity field along the line-of-sight should represent 
a signature of infalling material and may thus serve as 
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a rough guide to where shocks occur. In the left panel 
of Figure 9, we show the ratio of the gas temperatures 
for the full hydrodynamical simulation and the HPM 
simulation as a fuction of the velocity gradient of the 
gas. We first average the temperature in pixels within 
100 km/s from a minimum in the gradient of the pecu- 
liar velocity field in real space. Then, we average over 
the corresponding flux values, in redshift space. Before 
selecting the negative gradients in the hydrodynamical 
simulations, we have explicitely checked that the pecu- 
liar velocity fields are very similar in both simulations. 

The contour plots indicate the number density of 
points, which varies by a factor of 10 between adjacent 
contour levels. The bulk of the pixels in this panel is 
in regions with 5v/Sx ~ —0.5 and at 'hydrodynamical' 
temperatures that are about 10% lower than the cor- 
responding temperatures of the HPM simulation. The 
simulation at z = 2 shows a significantly increased 
amount of pixels at 5v/5x ~ — 1 with HPM tempera- 
tures that are colder than the corresponding tempera- 
tures in the hydrodynamical simulation. These differ- 
ences in the temperatures have an important effect on 
the simulated flux. 

In the right panel of Figure 9, we plot the differences 
in the flux for the same regions of infalling gas. Since 
the flux is observed in redshift space we have averaged 
the flux within 100 km/s velocity bins. There are no ob- 
vious trends for smaller values of Sv/Sx (i.e. "stronger" 
shocks). This is due to the fact that stronger shocks 
have a more complex temperature and density struc- 
ture which in the hydro simulation is represented more 
faithfully than in the HPM simulation. As a result, the 
differences in the temperatures and fluxes actually tend 
to be averaged out for strong shocks. The scatter for 
positive values of the gradient Sv/6x also shows smaller 
scatter, both in the temperatures ratio and in the flux 
differences. Most of the differences at z = 2 arise from 
regions of infalling gas that are not modelled accurately 
by the HPM method. This suggests that at least part of 
the discrepancy at low redshift is due to the increased 
amount of shock-heated gas probed by the Ly-a forest 
at lower redshift. 

3.3 The effect of the numerical parameters H B 
and A s 

As discussed in section 2.2, we need to specify the pa- 
rameters H s and A B which describe the smoothing of 
the gas density and of the gravitational force field in 
the HPM simulations. There is no obvious optimum 
choice for these parameters, so a choice needs to be 
made by comparing to the full hydrodynamical simu- 
lations. For changes of H s , which controls the smooth- 
ing of the pressure field, the resulting differences are 
very small at large scales (less than 1%). They are only 
weakly scale- and redshift-dependent and only slightly 
increase at small scales for H B in the range 1.5-3. We 
have therefore fixed H s = 3 for all simulations. 

Varying the parameter A s , which controls the 
smoothing of the gravitational force field, has a some- 
what larger effect. In Figure 10, we show the differences 
between the flux power spectrum of a HPM (JV gr id = 



600) simulation and that of the full hydrodynamical 
simulation for different values of A s , at three different 
redshifts z = 2,3, 4. The differences are typically a few 
percent but can be as large as 10%. It is not obvious 
which value for A s represents an optimum choice. One 
possibility is to impose a certain requirement for the 
maximum allowed force anisotropy generated by a point 
mass in the PM scheme. Using such a criterion, we have 
set A B — 1.25, which gives typical PM force errors less 
than 1 per cent. 

3.4 CPU time and memory requirements 

In Table 1, we summarise the total CPU time (in wall- 
clock seconds) required by the simulations to run to 
z — 2, and their memory requirement (in Gbytes). We 
include simulations with a range of particle numbers 
and resolutions of the PM mesh, all for a box size of 
30 Mpc/ft. The HPM simulation with jV gI .i d =600 has 
run about 20 times faster than the hydrodynamical 
SPH/TreePM simulation at the corresponding resolu- 
tion, but has a three times larger memory requirement. 
The A?g r id = 1200 HPM simulation, which as we saw gave 
a good agreement with the full hydrodynamical simula- 
tions in terms of the gas and dark matter distribution, 
is still faster than the SPH simulation by a factor 10, 
but its memory requirement is very large. We note that 
the simulations with a very high resolution of the PM 
mesh (iVg r id = 1200) have been difficult to run because 
of their very large memory requirement, which is close to 
the total amount available on the COSMOS computer 
we used. 



4 DISCUSSION AND CONCLUSIONS 

We have compared full hydrodynamical simulations car- 
ried out with the SPH/TreePM code GADGET- 2 to sim- 
ulations that used the HPM method. The latter scheme 
was implemented by us in GADGET-2 , and we com- 
pared this implementation with the independent code 
of GH. Our comparison was performed at redshifts 
2 = 2, z = 3 and 2 = 4. Our main results can be sum- 
marised as follows. 

• The dark matter and gas distributions of HPM sim- 
ulations with GADGET-2 converge well to the full hy- 
drodynamical simulations with the SPH/TreePM code. 
For a PM mesh with > 6 3 more mesh cells than par- 
ticles in the SPH simulations, the difference in the pdf 
of the gas and matter distributions are less than 1 per- 
cent. The same is true for the matter power spectrum 
at wavenumbers up to 20 times the fundamental mode 
of the box for a mesh with 1200 3 . At smaller scales the 
differences in the power spectra strongly increase due to 
lack of resolution of the HPM grid. 

• The pdf of the flux distribution of HPM simulations 
with GADGET-2 does not converge to that of the full 
hydro simulations. At low levels of absorption (F > 0.8) 
the differences (10% and more) do not decrease with 
increasing resolution at all three redshifts examined. At 
2 = 2, there is an additional large difference at low 
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Figure 8. Fractional differences between: 1) an HPM simulation with 2 X 200 3 and with 7V grid = 1200 and a full SPH 
hydrodynamical simulation (dotted line); 2) between a SPH simulation with zero artificial viscosity and with a superimposed 
temperature-density relation of the HPM runs and a full SPH hydrodynamical simulation (dashed line). Results are shown at 
z = 2,3,4 in the left, middle and right panels, respectively. Spectra have been scaled to reproduce the same effective optical 
depth. 
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Figure 9. Role of shock-heated gas. Left: Ratio of the temperatures in simulations run with the full hydro and HPM imple- 
mentation with A r gr ; t j=600 with the same initial conditions, plotted as a function of the gradient of the peculiar velocity field 
along the line of sight. Right: Differences in the simulated flux values. The contour plots represent the number density of points 
in the 2D plane and the number density increases by an order of magnitude at each contour level. The dashed, continuous and 
dotted lines are for z = 2, 3,4, respectively. 



flux levels which rises to 50% at the lowest flux levels. 
The latter difference is most likely due to the larger 
proportion of absorption by dense shock-heated gas at 
z — 2 which is not modelled well by the HPM method. 

• At redshifts z = 3 and z — 4, the flux power spec- 
trum of HPM simulations with GADGET- 2 does con- 
verge to that of the full hydrodynamical simulations 
up to a scale independent offset. For a HPM simula- 



tion with box size of 30/i 1 Mpc and a PM mesh with 
1200 3 cells this offset is about 5% -7% at wave numbers 
0.002 s/km < k < 0.05 s/km. At z = 2, however, there 
are large scale-dependent differences between the flux 
power spectrum of the HPM simulation and the full hy- 
drodynamical simulation which are as large as 20-40%. 
These differences are again most likely due to the larger 
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Figure 10. Effect of the smoothing parameter A B of the HPM implementation in GADGET-2 for simulations of a 30 Mpc//i 
box with 2 X 200'^ particles. We plot the fractional differences between the ID flux power spectra of a model with A B set to 1.25 
(thin line with filled triangles), 1.1 (thin line with empty triangles), 1 (thick line), 0.9 (thin line) and the full hydrodynamical 
simulation. The results are shown for three different redshifts z = 2, 3, 4 as dashed, continuous and dotted curve, respectively. 
Note that the flux power spectra have not been scaled to reproduce the same effective optical depth. 



proportion of absorption by dense shock-heated gas at 

2 = 2. 

• The HPM implementation of GADGET-2 and the 
code by GH give similar results (to within a few percent) 
for the same initial conditions, provided a slightly higher 
resolution of the PM grid is used for GADGET-2 . This 
offset is a result of the PM force smoothing done by 
GADGET-2 , which is adjustable. The results obtained 
above should thus hold in a similar form for the GH 
code. 

The HPM method involves two main simplifica- 
tions compared to full hydro- simulations, calculating 
the pressure in an approximate way and estimating the 
temperatures based on the density alone. The HPM 
approximation does a good job in modelling the gas 
and matter distribution on the scales relevant for the 
Lyman-a forest suggesting that the first approximation 
works well. The situation for an accurate prediction 
of the flux distribution is quite different and we have 
shown that the treatment of the thermal state in the 
HPM approximation is the main problem for accurate 
predictions of the flux distribution. The strong depen- 
dence of the transmitted flux on the thermal state of 
the gas together with the crude approximation of the 
thermal state in the HPM approximation leads to large 
and not always intuitive scale- and redshift-dependent 
differences in the flux distribution between HPM and 
the full hydrodynamical simulations. 

For the flux power spectrum, these differences are 
less important than for the pdf of the flux distribution. 
Our results suggest that at z = 3 and z — 4 the gain 
in speed offered by HPM simulations may still make 



them an attractive tool to obtain predictions of the flux 
power spectrum for a wide range of parameters. This 
will, however, require very careful calibration with full 
hydrodynamical simulations, and it appears doubtful 
that HPM simulations are suitable to model the depen- 
dence of the flux power spectrum on the thermal state 
of the gas accurately. The rather large memory require- 
ment of HPM simulations with sufficient resolution to 
reach convergence also partially offsets the advantage 
of their higher speed. Our results further suggest that 
at lower redshift the larger proportion of absorption by 
dense shock-heated gas makes HPM simulations unsuit- 
able for accurate predictions of the flux power spectrum. 

Currently the observational uncertainties regard- 
ing the thermal state of the IGM are still rather large. 
The results of quantitative studies of the matter power 
spectrum with Lyman-alpha forest data are therefore 
generally marginalized over a wide range of simple 
temperature-density relations. The difficulties of simple 
HPM implementations with modeling the effect of the 
thermal state accurately may therefore be less impor- 
tant than suggested by our discucussion so far. However, 
improved measurements of the thermal state of the gas 
utilizing the Doppler parameter distribution, the flux 
PDF and the small scale flux power spectrum are an 
important prerequisite for reducing the errors of mea- 
surements of the matter power spectrum from Lyman- 
a forest data. 

Accurate modeling of the thermal state of the gas 
will be required to take full advantage of an reduced un- 
certainty regarding the thermal state of the IGM. For 
HPM simulations this will almost certainly require a 
signifificant improvement of the modelling of the ther- 
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mal state, e.g. by introducing some scatter in the tem- 
perature density relation. Full hydrodynamical simu- 
lations could thereby be used to quantify and cali- 
brate this scatter and to investigate possible correla- 
tions of the scatter with physical quantities. Such mod- 
elling would obviously greatly benefit from more pre- 
cise observational estimates of the parameters describ- 
ing the temperature-density relation which may be pos- 
sible with the use of the flux power at smaller scales 
and from an estimate of the scatter in the temperature 
density relation using higher order statistics such as the 
bispectrum (Mandelbaum et al. 2003, Viel et al. 2004, 
Fang & White 2004). It will then also be important (in 
HPM and full hydro simulations) to model other physi- 
cal aspects affecting the thermal state of the gas as the 
presence of galactic winds and temperature/UV fluctu- 
ations due to the reionization of Hell. 
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